#########################
### Results for Justice as Checks and Balances
### Using Panel Data 
##########################

### This script works with panel data at province-decade level
##  This code creates:
##  Determinants for number of files (First stage bias): Figure 4 and Tables A4.1 and A4.3
##  Favorable Ruling and Welfare Effects (S2)
##  Claims and rebellion (S6)

# Preliminaries --
rm(list=ls())
gc()

### Libraries
###########################
# Relevant libraries
# Next, we download packages that H2O depends on.
pkgs <- c("plm", "dplyr", "magrittr", "ggplot2", "lmtest", "stargazer",
          "aod", "multiwayvcov", "mediation", "DataCombine", "graphics", "lmtest", "xtable"
)
for (pkg in pkgs) {
  if (! (pkg %in% rownames(installed.packages()))) { install.packages(pkg) }
  lapply(pkg, library, character.only=T )
}


#### Read panel data 
agn_panel_pob_decades <- read.csv("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/agn_panel_2021.csv")

#### merge with high merit cases 
high_merit <- read.csv("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/high_merit_panel.csv") %>% dplyr::select(name, decade, high_merit)
high_merit <- high_merit %>% dplyr::select(name, decade, high_merit)
agn_panel_pob_decades <- agn_panel_pob_decades %>% left_join(high_merit) %>% mutate(high_merit_log=log(high_merit),
                                                                                    high_merit_per=high_merit/files_ind,
                                                                                    high_merit_per=ifelse(is.na(high_merit),0,high_merit_per)
                                                                                    )                                                                                    
### Read rebellions data 
rebel <- read.csv("C://Users/franc/Dropbox/PhD/dissertation/strategies_resistance_accomodation/rebel_2021.csv")
rebel <- rebel %>% dplyr::select(-X)
### Merge with panel 
agn_panel_pob_decades <- agn_panel_pob_decades %>% left_join(rebel)

#######################
##  
##   Structural bias analysis
##   Table A4.1, Table 4.3, and Figure 4
#######################

##### Model by decades
model2_files_d <- lm(log_files  ~ pob_ln + pop_decrease + power + cat2+ Fragmented + franc+ multiling +
                       distance+small_scale2+ as.factor(decade), data=agn_panel_pob_decades)  

#summary(model2_files_d)
### QuasiPoisson
model2_files_d_qp <- glm(files_ind  ~ pob_ln + pop_decrease +power + cat2+ Fragmented + franc+ multiling +
                           distance + small_scale2+ as.factor(decade), data=agn_panel_pob_decades, family="quasipoisson")

####Point estimates plot
model2_files_d_coef  <- coeftest(model2_files_d )
model2_files_d_coef  <- model2_files_d_coef [c(2:10), c(1:2) ]
model2_files_d_coef  <- as.data.frame(model2_files_d_coef )
colnames(model2_files_d_coef ) <- c("Est", "SE")
model2_files_d_coef $name <- c("Indigenous pop. (log)", "Indigenous pop. declining",
                               "Strong local elites", "Extractive potential", "Fragmented province", 
                               "Franciscan jurisdiction", "Multilingual", "Distance to Mexico City", "Protests")
#### Order
model2_files_d_coef$ord <- c(1:nrow(model2_files_d_coef ))

###########################
##  Figure 4
## Coefplot
##########################
g <- ggplot(model2_files_d_coef , aes(x=reorder(name, -ord), y=Est, ymin=Est-(1.96*SE), ymax=Est+(1.96*SE))) + geom_pointrange()
g <- g + geom_hline(yintercept=0, color="black", linetype="dashed")
g <- g + theme_bw()+ coord_flip() + ylab("Coefficient") +xlab("")
g <- g + ggtitle("") + theme(axis.text=element_text(size=12),
                                                                 axis.title=element_text(size=12,face="bold"))
g
### Save 
ggsave("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/FINAL_version_2021/figures_paper/figure4.tiff",g, width=5)

##########
#### Table A4.3
#### Only high merit cases 
##### Model by decades
model2_files_high <- lm(high_merit_per ~ pob_ln +pop_decrease +power + cat2+ Fragmented + franc+ multiling +
                       distance +small_scale2 + as.factor(decade), data=agn_panel_pob_decades)  

model2_files_high_qp <- glm(high_merit_per  ~ pob_ln + pop_decrease +power + cat2+ Fragmented + franc+ multiling +
                           distance + small_scale2+ as.factor(decade), data=agn_panel_pob_decades, family="quasipoisson")


## Table A4.1
stargazer(model2_files_d ,model2_files_d_qp ,
          #column.labels = c("All Cases", "High-Complexity Cases"),
          #column.separate = c(2, 2),
          dep.var.labels= c("OLS", "QP"),
          covariate.labels =c("Indian pop. (log)", "Pop. Decreasing", "Strong LE",
                              "Category", "Fragmented", "Franciscans", "Multilingual",
                              "Distance to Mexico City" , "Protests"),
          keep = c("pob_ln"  ,"pop_decrease", "power"  ,"cat2" ,"Fragmented" , "franc","multiling" ,
          "distance" , "small_scale2"), 
          keep.stat = c("n", "rsq"),
          #out="bias_table.tex", 
          label="bias", 
          add.lines=list(c("Time FE", "Y", "Y", "Y", "Y")),
          title = "Determinants of Claims")

## Table A4.3
## Table high merit
stargazer(model2_files_high ,model2_files_high_qp  ,
          #column.labels = c("All Cases", "High-Complexity Cases"),
          #column.separate = c(2, 2),
          dep.var.labels= c("OLS", "QP"),
          covariate.labels =c("Indian pop. (log)", "Pop. Decreasing", "Strong LE",
                              "Category", "Fragmented", "Franciscans", "Multilingual",
                              "Distance to Mexico City" , "Protests"),
          keep = c("pob_ln"  ,"pop_decrease", "power"  ,"cat2" ,"Fragmented" , "franc","multiling" ,
                   "distance" , "small_scale2"), 
          keep.stat = c("n", "rsq"),
          #out="bias_table.tex", 
          label="bias", 
          add.lines=list(c("Time FE", "Y", "Y", "Y", "Y")),
          title = "Determinants of Claims (High merit cases)")



#######################
##  
##   Welfare Effects (supplementary material S2)
##   Effect of protection on population decline
##   Table S2.1
#######################

#### Effect of protection on:
#### Population decrease =1
#### Current decade
model_protect <- lm(pop_decrease ~ protect_per +pob_ln + as.factor(decade) + name, 
                    data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]
)
model_protect_vcov <- cluster.vcov(model_protect , ~ agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]$name)
model_protect_vcov <- coeftest(model_protect , model_protect_vcov)

### past decade (lag)
model_protect_1 <- lm(pop_decrease ~ protect_per_1 +pob_ln+ as.factor(decade)+name, 
                      data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),])
model_protect_1_vcov <- cluster.vcov(model_protect_1  , ~ agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]$name)
model_protect_1_vcov <- coeftest(model_protect_1  , model_protect_1_vcov )

### Acumulated (present and past decade)
model_protect_acum <- lm(pop_decrease ~ protect_per + protect_per_1+pob_ln+ as.factor(decade)+name, 
                         data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),])
model_protect_acum_vcov <- cluster.vcov(model_protect_acum , ~ agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]$name)
model_protect_acum_vcov <- coeftest(model_protect_acum , model_protect_acum_vcov )

#### First differences
### Change totals
model_protect_first <- lm(pop_change ~  protect_change+as.factor(decade)+name, 
                          data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),])

model_protect_first_vcov <- cluster.vcov(model_protect_first, ~ agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]$name)
model_protect_first_vcov <- coeftest(model_protect_first , model_protect_first_vcov)


### Percentages change
model_protect_first_per <- lm(pop_change ~  protect_per_change+as.factor(decade)+name, 
                              data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),])

model_protect_first_per_vcov <- cluster.vcov(model_protect_first_per, ~ agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"),]$name)
model_protect_first_per_vcov <- coeftest(model_protect_first_per , model_protect_first_per_vcov)

##### units
units <- length(unique(agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico")]))

##### Table S2.1
#### Favorable ruling and population change
stargazer(model_protect , model_protect_1,model_protect_acum, model_protect_first ,model_protect_first_per ,
          se   = list(model_protect_vcov[,2] , model_protect_1_vcov[,2],model_protect_acum_vcov[,2], 
                      model_protect_first_vcov[,2], model_protect_first_per_vcov[,2]
          ) ,
          p    = list(model_protect_vcov[,4] , model_protect_1_vcov[,4],model_protect_acum_vcov[,4],
                      model_protect_first_vcov[,4], model_protect_first_per_vcov[,4]
          ),
          
          dep.var.labels= c("Population Decrease",  "Population Change"),
          covariate.labels =c("Favorable Ruling",  "Favorable Ruling t-1",  "Favorable Ruling Change","Favorable Ruling Percentage Change"
          ),
          #column.labels = c("Pre Bourbon Reform", "Post Bourbon Reform"),
          #column.separate = c(1,1),
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("protect", "protect_1", "protect_change", "protect_per_change"), 
          keep.stat = c("n", "adj.rsq", "ll"),
          add.lines=list(
            c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y"),
            c("Time FE", "Y", "Y", "Y", "Y", "Y", "Y"),
            c("Units", rep(units, 6))
          ),
          label="welfare_effects",
          title =("Favorable Ruling and Population Change") 
          #out="welfare.tex"
          )


#######################
##  
##   Rebellion Analysis (supplementary material S6)
##   Tables S6.2 and S6.3
#######################

#####  Total large scale 
model_reb1      <- lm(rebel_hist2 ~ protect+ files_ind  + pob_ln + as.factor(name) + as.factor(decade), data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico") ,] )
model_reb1_vcov <- cluster.vcov(model_reb1 , ~ agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico") ])  ## clustered se
model_reb1_vcov <- coeftest(model_reb1  , model_reb1_vcov)

### Total Small scale events 
model_reb2  <- lm(small_scale2~ protect+ files_ind +  pob_ln + as.factor(name) + as.factor(decade) , data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico"), ]) 
model_reb2_vcov <- cluster.vcov(model_reb2 , ~ agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico") ])  ## clustered se
model_reb2_vcov <- coeftest(model_reb2  , model_reb2_vcov)


#####  Total large scale 
model_reb3      <- lm(rebel_hist2 ~ protect_per  + pob_ln + as.factor(name) + as.factor(decade), data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico") ,] )
model_reb3_vcov <- cluster.vcov(model_reb3 , ~ agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico") ])  ## clustered se
model_reb3_vcov <- coeftest(model_reb3  , model_reb1_vcov)


### Total Small scale events 
model_reb4  <- lm(small_scale2~ protect_per +  pob_ln + as.factor(name) + as.factor(decade) , data=agn_panel_pob_decades[which(agn_panel_pob_decades$name!="Mexico") ,]) 
model_reb4_vcov <- cluster.vcov(model_reb4 , ~ agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico") ])  ## clustered se
model_reb4_vcov <- coeftest(model_reb4  , model_reb4_vcov)


##
units <- length(unique(agn_panel_pob_decades$name[which(agn_panel_pob_decades$name!="Mexico") ]))

######## Stargazer
##### Table A12
stargazer(model_reb1,model_reb2, model_reb3, model_reb4,      
          se   = list(model_reb1_vcov[,2] , model_reb2_vcov[,2],model_reb3_vcov[,2], 
                      model_reb4_vcov[,2]
          ) ,
          p    = list(model_reb1_vcov[,4] , model_reb2_vcov[,4],model_reb3_vcov[,4],
                      model_reb4_vcov[,4]
          ),
          dep.var.labels= c("Rebellions", "Riots, Revolts, etc.","Rebellions", "Riots, Revolts, etc."),
          covariate.labels =c("Successful Claims (tot.)",  "% Successful Claims"
          ),
          #column.labels = c("Total Susccesful Claims", "Proportion of Succesful Claims"),
          #column.separate = c(2,2),
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("protect", "protect_per"), 
          keep.stat = c("n", "adj.rsq", "ll"),
          add.lines=list(
            c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y"),
            c("Time FE", "Y", "Y", "Y", "Y", "Y", "Y")
            #c("Units", rep(units, 6))
          ),
          label="rebel_effects",
          title =("Uprisings and Successful Claims") 
          #out="rebel.tex"
          )

#### Table S6.3
### causality test between rebellion and cases
### Does rebellion causes claims?
r_to_c <- grangertest(agn_panel_pob_decades$files_ind~ agn_panel_pob_decades$small_scale2, order=1) # not significant
### Does claims cause rebellion?
c_to_r <- grangertest(agn_panel_pob_decades$small_scale2~ agn_panel_pob_decades$files_ind, order=1) ## More likely

### Does rebellion causes protection?
r_to_p <- grangertest(agn_panel_pob_decades$protect_per~ agn_panel_pob_decades$small_scale2, order=1)
### Does protection cases rebellion?
p_to_r <- grangertest(agn_panel_pob_decades$small_scale2~ agn_panel_pob_decades$protect_per, order=1)

#### Table
y <- c("Total files", "Small scale events", "% protection", "Small scale events")
x <- c("Small scale events", "Total files", "Small scale events", "% protection")

p_val <- c(r_to_c[2,4],c_to_r[2,4], r_to_p[2,4], p_to_r[2,4])

tab <- cbind.data.frame(y,x,p_val)
print(xtable(tab),include.rownames=FALSE)


